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Department  of  Mathematics,  University  of  South  Carolina,  Columbia,  SC  29208 


Abstract —  Two  elementary  algorithms  are  introduced  for  image  compression  each  of  which 
is  based  on  efficient,  lossless  encoding  of  quantized  bi-orthogonal  wavelet  coefficients.  Application 
of  this  type  of  algorithm  is  applied  to  several  standard  test  images  using  regular  and  hyperbolic 
wavelet  bases,  and  comparisons  are  given  to  Shapiro’s  EZW  algorithm.  Peak  signal  to  noise  ratio 
improvements  typically  of  0.6— 0.8  dB  are  demonstrated.  Generalizations  of  this  type  of  algorithm  to 
non-square  images  and  higher  dimensions  are  also  briefly  described. 


1.  Introduction 

Image  compression  is  important  in  many  applications,  such  as  image  transmission,  feature 
extraction  and  denoising  of  data  [6,  8,  11,  3].  Our  original  intent  was  to  develop  reliable  data 
compression  algorithms  for  three-dimensional  time-dependent  scalar  and  vector  fields  with  gen¬ 
eral  index  limits  (i.e.  non-dyadic).  In  particular,  a  requirement  for  3D  data  compression  arose 
in  our  development  of  utilities  for  the  interactive  tracking  and  steering  of  remote  simulations  over 
small  bandwidth  or  congested  communication  connections  [14,  15].  The  underlying  rationale  is 
that  scientists  and  practitioners  need  the  capability  to  quickly  assess  a  simulation  and  correspond¬ 
ingly  terminate  or  alter  its  progress.  Interactive  work  on  large  scale  simulators  typically  require 
data  compression  algorithms  which  permit  a  careful  analysis  of  three  dimensional  data  at  various 
frequencies  and  resolutions.  Wavelet  analysis  is  a  natural  framework  for  these  studies. 

In  our  investigations,  we  developed  various  encoding  algorithms  for  multi-dimensional  data  and 
tested  them  in  the  two-dimensional  case  by  comparing  the  performance  with  existing  algorithms 
for  2-D  images.  We  obtained  the  somewhat  surprising  results  that  these  algorithms  performed  very 
well  in  the  tests,  although  they  were  much  more  simple  to  implement  and  generalize.  The  basic 
steps  of  our  image  compression  algorithm  include  (1)  application  of  discrete  bi-orthogonal  wavelet 
transformations  to  an  extension  of  the  image  data,  (2)  a  one-step  quantization  according  to  the 
precision  required,  (3)  a  natural  ordering  of  the  coefficients  by  index,  (4)  a  preprocessing  procedure 
of  encoding  the  space-frequency  correlations  of  the  coefficients,  followed  by  (5)  an  application  of 
a  Q-coder  algorithm.  The  reconstruction  (to  within  quantization)  involves  the  simple  inversion  of 
these  operations  in  the  reverse  order. 

Steps  1  and  3-5  are  generally  lossless,  but  in  practice  noise  may  be  present  in  the  original  image. 
The  information  loss  occurs  in  the  quantization  stage,  which  may  be  estimated  by  approximation 
results  of  DeVore  et  al.  (see  [7]  and  inequality  (21)  below).  As  we  show,  the  combined  algorithm 
is  very  easy  to  implement  and,  compared  to  other  algorithms,  is  computationally  efficient  even  for 
multidimensional  data. 

1This  work  was  supported  in  part  by  Martin  Marietta  Subcontract  SK966V,  DOE  Grant  DE-FG05-95ER25266, 
ONR  Grant  N00014-96- 1-1003,  and  DoD  Grant  N00014-97- 1-0806. 
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In  the  following  sections  we  will  discuss  each  of  these  steps.  Section  2  outlines  the  multiresolu¬ 
tion  analysis  for  images  based  on  tensor-products.  Section  3  discusses  the  quantization  and  error 
estimate,  while  Section  4  describes  the  ordering  of  the  quantized  data  without  storing  the  location 
of  each  datum.  We  have  developed  preprocessing  algoritms  which  have  many  variants,  two  of  these 
are  described  in  Section  5.  In  Section  6,  we  present  the  results  of  computational  experiments  with 
these  encoders  and  compare  their  effectiveness  with  published  results  of  Shapiro  and  an  implemen¬ 
tation  of  that  algorithm.  We  provide  concluding  remarks  and  observations  in  Section  7  which  relate 
the  encoding  technique  with  entropy  and  data  complexity.  A  preliminary  version  of  this  research 
was  reported  earlier  in  [12].  We  would  like  to  express  our  thanks  to  Ronald  DeVore  for  discussions 
concerning  the  results  of  this  paper. 

2.  Wavelet  Transforms 

In  this  section  we  briefly  describe  biorthogonal  locally  supported  wavelets  and  the  corresponding 
multiresolution  analysis  (see  [2,  16,  22]  for  details  and  additional  references).  We  begin  with  an 
outline  of  the  development  in  the  univariate  case  with  an  emphasis  on  how  to  implement  both 
decomposition  and  reconstruction  of  discrete  data,  and  complete  the  section  by  describing  2D 
tensor  product  bases  for  images. 

2.1  Multiresolution  analysis. 

A  multiresolution  analysis  of  L2(R)  is  defined  as  a  ‘ladder’  of  closed  subspaces  14  which  satisfy 
the  following  properties: 

(l.a)  There  exists  a  scaling  function  p  with  a  non- vanishing  integral,  so  that  the  closure  of  the 
span  of  its  integer  translates  ( Pj(x )  :=  <p(x  —  j))  is  Vo-  For  computational  purposes  we 
will  assume  that  p  has  compact  support.  Moreover,  we  assume  the  collection  {pj \j  £  Z} 
is  a  Riesz  basis  for  Vq:  there  are  constants  C\ .  C2  so  that  if  /  £  Vo,  then  /  = 
and 

^l||fe}||^<||/||L2(R)<C'2||{c,-}||,2 


(1.6) 

Vk  c  14+1, 

keZ 

(l.c) 

(J  Vfc  =  L2(R)  , 

Z) 

11 

fee  z 

fcez 

(1  -d) 

f{x)  £  Vk  /( 2x)  £  14+i  for  all  k  £  Z 

Although  some  of  these  properties  are  redundant,  this  formulation  is  best  for  our  purposes.  For  a 
given  scaling  function  p,  property  (l.b)  (with  k  =  0)  implies  that  a  refinement  relation  holds 

(2)  p(x)  =  2Y^otjp(2x  -  j). 

j 
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By  property  (l.a)  the  coefficients  must  be  unique.  We  assume  further  that  the  coefficients  ay  satisfy 
the  property 

ai  =  1  • 
iez 

In  practice,  one  usually  determines  coefficients  with  desirable  properties  for  the  resulting  scale  (and 
wavelet)  functions  and  show  that  the  above  properties  are  satisfied.  As  contrasted  with  the  usual 
orthogonal  development  of  wavelets,  the  theory  of  biorthogonal  wavelets  instead  makes  use  of  two 
ladders  of  subspaces,  one  for  the  primal  space  and  another  for  its  dual.  In  the  case  of  Hilbert  spaces 
L2,  these  spaces  may  be  regarded  as  coinciding. 


2.2  Univariate  bi-orthogonal  decompositions. 

A  biorthogonal  wavelet  analysis  begins  with  two  dual  scale  functions  p,  p  which  each  provide  a 
multiresolution  analysis  14, 14 >  respectively  of  L2 .  Actually,  one  should  think  of  Vy.  as  providing  a 
decomposition  of  L2  and  Vy-  performing  a  similar  role  for  its  dual.  We  denote  by  ay  the  coefficients 
following  from  property  (2)  for  the  scale  function  p.  We  define  corresponding  dual  wavelet  functions 
by 

(3. a)  ip(x)  =  2^2Pjip(2x  -  j)  ,  ^(x)  =2j2Pj‘P(2x  -  j) 

j  j 

where 

(3.6)  (3j  =  (-l)Jdi _ y  ,  fa  =  (-l)Jai_y. 

We  also  assume  that  the  dual  scale  functions  p,  p  have  the  following  properties: 


(4) 


then  it  follows  that 


lR 


p(x)  ip(x  —  j)  dx 


$o j 


(5) 


t  R 


iAx )  dx 


where  we  are  using  the  standard  notation 

(6)  4>k,j(x)  =  V¥(, fi(2kx  -  j ) 


for  the  L2-normalized  translated-dilates  of  a  function  (f>.  Here  the  index  k  denotes  the  scale  and  j 
indicates  the  integer  shift.  For  our  computational  purposes  we  assume  that  each  of  the  refinement 
sums  in  (2)  is  finite  with  an  even  number  of  symmetric  ‘filter’  coefficients  ay,  ay: 


(7)  P{x)  =  2  J2  ajP(2x~j),  <p(x)  =  2  ajp{2x-j). 

j  =  l—U  j  =  l—u 

The  scale  and  wavelet  functions  p  and  if  generate  corresponding  spaces  14,  IVy.  defined  by  the  closed 
linear  spans 

Vk  =  spanjy^y  :  j  G  Z} 


3 


Wk  =  spanjV’fcj  :j€  Z} 

Similarly,  we  denote  by  14  ,Wk  the  collection  of  spaces  corresponding  to  (p  and  4,  respectively.  It 
follows  from  the  properties  of  the  dual  scale  functions  and  the  definitions  of  the  dual  wavelets  that 

(8)  I4+i  =  Vk  ©  Wk,  Vk+i  =  Vfc  ©  Wfc. 

Although  Wk  may  not  be  the  orthogonal  complement  of  14,  the  following  relationships  hold: 

(9)  Wk  _L  Vk,  Vk  _!_  Wfc. 

The  decomposition  (8)  is  repeated  recursively  to  show  for  each  nonnegative  n  €  Z 

(10)  14  =  I4-i  ©  W„_i 

=  14-2  ©  (wn_2  ©  wn_i) 


=  Vo  ©  (w0  ©  •  ■  •  ©  wn_i)  . 


Each  /  G  1 4  may  then  be  written  uniquely  as 


(11) 


/ 


cj  Pnj 

ie  Z 


c7  Pn-lj  +  V’n-lJ 

iez  jez 


2<pn-2,j  + 

3  ez 


24n-2j 

J'ez 


+  Vn-lj 

J6Z 


n—  1 


=  E  ci  Wj  +  E  ^ 


iez 


fc=o  \jez 


which  is  called  its  biorthogonal  wavelet  decomposition. 


2.3  Decomposition  stage  for  discrete  functions. 

For  an  discrete  function  c  =  {c™ jjlo1  we  associate  the  function  /  £  14  defined  by 

/  =  £  C>nJ 

jez 


where  the  coefficients  are  extended  symmetrically  at  the  two  boundaries  j  =  0  and  j  =  2n  —  1  by 
the  formula 


(12) 


c-i-j,  3  =  -1,-2,  •••  ,-u 

c^+i-^p  j  =  2",  2"  +  1, . . . ,  2"  +  fi  -  1 

0,  otherwise 


in  order  to  allow  a  systematic  procedure  for  producing  the  next  coarser  level  scale  and  wavelet  co¬ 
efficients  as  one  procedes  recursively  through  each  step  represented  by  the  individual  lines  in  (11). 
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From  the  biorthogonality,  we  may  reproduce  c  through  ‘sampling’  /  by  integration  against  the  cor¬ 
responding  i finj.  In  particular,  as  k  =  n,  n  —  1, . . . ,  1  the  wavelet  and  scale  coefficients,  respectively, 
at  the  next  coarser  level  k  —  1  are  generated  by 

U  U 

(i3)  d) r1  =  s  E  a  4+i'  A1  =  ^  E  “i  4+i 

Z=l— u  1=1— u 


where  0  <  j  <  2k  1  —  1.  Similar  to  (12)  above,  we  use  the  symmetric  extension  of  the  coefficients 
at  each  level: 

j  =  -1,-2,..., -ft 
j  =  2k,2k  +  l,...,2k  +u-l 


(14) 


ckd  =  \  y-i’ 

c2k+1—j—V 
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to  provide  the  necessary  terms  so  that  the  process  may  proceed  level-by-level.  Symbolically  this 
may  be  represented  by  the  following  diagram: 


{c]-1}  — 

I 

1 

T 

\ 

\ 

K"1} 

K_1> 

2.4  Reconstruction  stage  for  discrete  functions. 

A  discrete  function  is  reconstructed  as  represented  by  the  following  diagram  by  moving  from 
the  coarse  to  the  fine  levels  (k  =  0, 1, . . ,  n  —  1): 


/ 

K°} 

using  the  reconstruction  filters 


rk  — 
c2j  - 


(15) 


-2j+l 


{A-1}  _>  {cn 


K"1} 


{cr1} 


K"  > 


=  \pl 


"kzi+  E  fa  dti 


=  \pl 


l  Lf  J  LSI 

E  c‘ 

V=-L^J  " 

( 

E  a2z+i  ckzi  +  e  fa+i  dj~i 

V=-Lfi 


Lnr-J 

E 

* — L§J 


Kl 


for  the  range  of  j  from  0  to  2k~l  —  1.  Here  again  we  take  for  the  discrete  function  {cfc_1}  its  even 
symmetric  extension  according  to  (14),  but  for  the  wavelet  coefficients  {d^-1}  we  must  use  the  odd 
symmetric  extension. 


2.5  Images  and  Tensor  Products  of  Wavelet  Bases 

A  square  image  can  be  represented  by  its  pixel  values  ctj,  ( i,j  =  0, 1, . . . ,  2n  —  1)  where  all  these 
values  are  integers  in  the  interval  [0,  255] .  If  we  write  these  pixel  values  in  a  matrix  Mq  and  apply 
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the  scaling  and  wavelet  transforms  (13)  to  each  row  of  Mo,  we  obtain  matrices  M\  and  M2  of  scale 
and  wavelet  coefficients  in  the  x- variable  with  y  held  fixed: 

AIq  Mi  M2 

In  each  of  Ml  and  M2  we  arrange  the  coefficients  first  according  to  scale  (finest  to  lowest)  and 
then  by  position  (left  to  right).  Applying  the  univariate  wavelet  transform  next  to  the  columns  of 
Mo,  Mi,  and  M2,  we  obtain  six  additional  matrices  M3, . . . ,  Mg: 

Mq  Al\  M2 

(16)  M3  M4  M5  , 

Mg  M’j  Mg 

where  Mi ,  M2  have  dimension  2n  x  (2n  —  1 ) ,  M3 ,  Mg  have  dimension  {2n  —  1 )  x  2n ,  and  the  remainder 
have  dimension  (2n  —  1)  x  {2n  —  1).  Any  2n  x  2n  linearly  independent  basis  elements  form  a  basis 
for  a  two  dimensional  wavelet  transform.  Although  it  is  an  abuse  of  terminology,  we  say  that  the 
coefficients  are  linearly  independent  when  the  corresponding  collections  of  scaling  and/or  wavelet 
functions  are  linearly  independent.  Two  bases  generally  used  are  the  hyperbolic  and  regular  bases. 
To  form  the  hyperbolic  decomposition,  one  uses  all  the  coefficients  in  the  matrix  Mg  (i.e,  tensor 
products  of  wavelets  from  all  scales),  the  last  row  in  M5  (all  wavelets  in  x  tensored  with  the  coarsest 
scale  function  in  y),  the  last  column  in  M7  (all  wavelets  in  y  tensored  with  the  coarsest  scale  function 
in  x)  and  the  coefficient  in  the  lower  right  corner  of  M4  (product  of  the  coarsest  scale  functions  in 
x  and  y) .  The  regular  basis  decomposition  is  the  more  standard  one  employed  in  image  processing 
and  uses  all  the  coefficients  corresponding  to  the  functions  <p(x)ijj(y),  il>(x)<p(y),  i/>( x )i/>(y)  scaled 
appropriately  with  square  support,  (i.e.,  corresponding  elements  in  M5,  M7  and  Mg),  together  with 
the  additional  coarsest  scale  coefficient  in  the  lower  right  corner  of  M4.  These  two  bases  will  be 
compared  in  our  computational  experiments  in  Section  6. 

3.  Quantization. 

After  an  application  of  wavelet  transforms,  we  obtain  an  indexed  collection  of  real- valued  coef¬ 
ficients.  In  the  bivariate  case,  if  N  =  2n  and  a  function  /  is  represented  as 

N—l 

f{x,y)=  ^  c™j  ip(2nx  —  i)<p(2ny  —  j) 
i,j= 0 

then,  for  convenience,  we  may  write  it  as 

N2- 1 

(17)  f{x,y)  =  ak<t>k{x,y) 

k= 0 

where  the  functions  </>*.  are  chosen  to  form  either  a  hyperbolic  or  regular  basis,  and  are  normalized  in 
the  Lebesgue  space  Lp  for  a  specific  p  in  the  interval  (0,  00).  In  the  seminal  paper  ([7],  see  also  [6]) 
on  the  application  of  nonlinear  approximation  to  wavelet  compression,  DeVore  et  al  established, 
among  many  other  results,  that  if  /  belongs  to  the  Besov  smoothness  space  ( Lq(I ))  and  is 
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represented  as  in  equation  (17)  and  if  /  has  a  similar  representation  with  coefficients  ak  which 
satisfy 

(18)  \ak  -  ak\  <  51/q, 
then 

(19)  ||/  —  I\\lp  =  0(5a/2). 

Here  a  is  a  measure  for  the  smoothness  of  /  (see  [7]  for  a  complete  explanation). 

Based  on  this  estimate,  uniformly  scaled  integers  will  be  used  to  approximate  the  real-valued 
coefficients  in  this  representation  according  to  a  given  precision  e.  Let  [a]  denote  the  closest  integer 
to  a  real  number  a.  For  any  e  >  0,  we  have: 

a  —  2e  [a/ 2e]  <  e. 

If  we  apply  this  with  a  =  ak,  and  define  the  approximation  /  by 

JV2  — 1 

(20)  f(x,y)  =  2e  [afc/2e]  tk(x,y) 

k=0 

then  the  coefficients  ak  of  /  satisfy  inequality  (18)  and  therefore  the  approximation  result  given 
in  (19)  implies  that 

(21)  \\f_f\\LP=0(e^/2). 

This  process  is  referred  to  as  quantization  and  e  is  called  the  threshold.  Note  that  all  coefficients  of 
/  whose  absolute  value  is  smaller  than  the  threshold  are  quantized  to  be  zeros. 

4.  Ordering  of  Coefficients. 

Many  compression  algorithms  are  based  upon  the  assumption  that  a  majority  of  the  data  after 
quantization  are  zero.  In  [7]  an  algorithm  was  described  which  discarded  quantized  coefficients 
whose  magnitude  was  smaller  than  the  threshold.  The  remaining  quantized  coefficients  were  ar¬ 
ranged  in  decreasing  order  and  their  indices  were  recorded.  A  careful  description  was  presented  of 
the  types  of  quantization  at  each  scale  level  and  methods  for  packing  the  location  of  the  quantized 
coefficients  to  reduce  the  number  of  bits  per  pixel  required  to  transmit  the  compressed  image. 

We  avoid  the  overhead  of  ordering  the  coefficients  and  attempt  to  minimize  the  effort  in  the 
transmission  of  the  locations  of  the  significant  coefficients.  Our  approach  uses  the  observation  that 
the  data  are  naturally  arranged  in  a  rough  increasing  order  when  viewed  relative  to  their  locations 
within  the  matrices  described  in  Section  2.4.  This  rough  estimate  for  a  given  function  /  is  indicated 
by  the  following  elementary  inequality 

(22)  f  ft  <\b-a\1/2\\f-P\\2 


7 


where  (j)  is  normalized  in  L2,  has  m  vanishing  moments  and  is  supported  in  [a,b\  and  P  is  any 
polynomial  of  degree  at  most  m  on  [a,  b\.  The  inequality  follows  from  the  standard  approximation 
argument.  Since  (j)  has  m  vanishing  moments,  /  P  cf>  =  0  and 


/(/  -  P)4> 

\b-a\l/q 


<11/- 
—  P II 

1  OO* 


P\\Li[a,b]  ll/ll Lp 


Roughly  speaking,  inequality  (22)  says  that  the  smaller  the  support  of  of  a  wavelet,  the  smaller 
the  coefficient  of  /  corresponding  to  that  wavelet.  In  the  matrix  of  coefficients  (16),  however,  the 
support  of  each  coefficient  is  directly  related  to  their  matrix  positions,  with  the  upper  left  corner 
coefficients  having  the  smallest  support  and  those  in  the  lower  right  corner  the  largest.  Therefore, 
as  we  index  through  the  appropriate  sub-matrices  (for  a  given  basis),  we  move  from  smaller  to 
larger  supports  to  ensure  that  generally  coefficients  with  small  values  occur  first. 


5.  Encoding  of  Quantized  Coefficients. 

In  this  section,  we  describe  two  relatively  simple  approaches  for  encoding  the  quantized  wavelet 
coefficients  resulting  from  the  decomposition  procedure  described  in  the  last  section.  We  use  to 
our  advantage  the  spatial-frequency  localization  of  features  to  further  condense  the  significant 
information  extracted  by  the  decomposition  procedure.  The  methods  we  describe  are  both  lossless 
encoders  which  pre-process  the  data  to  enhance  the  subsequent  application  of  standard  encoders, 
such  as  Q-coder  [17],  for  example.  The  first  approach,  which  we  refer  to  as  “interleaving,” 
compresses  the  rows  (or  columns)  of  the  array  of  wavelet  coefficients  by  recursively  interleaving 
adjacent  coefficients  in  each  row  of  the  coefficient  matrix  and  only  retaining  the  significant  bits. 
The  second  approach,  refered  to  as  “bit-stream  encoding,”  makes  use  of  the  higher  correlations 
which  exist  in  multivariate,  logically  rectangular  data. 

5.1  The  Interleaving  Method. 

Many  wavelet-based  compression  algorithms  rely  upon  the  assumption  that  a  majority  of  the 
quantized  data  vanish.  Improving  the  approximation  to  /  requires  decreasing  the  threshold  to 
produce  additional  nonzero  coefficients.  Nevertheless,  many  of  the  quantized  coefficients  remain 
zero,  and  the  ones  no  longer  zero  will  be  small  integers. 

The  first  encoding  method  uses  this  empirical  observation  and  the  fact  that  small  signed  integers 
require  a  small  number  of  bits  for  their  representation.  The  algorithm  condenses  pairs  of  these 
integers  to  encoded  integers,  using  the  natural  ordering  within  the  matrices  (16)  as  described  in 
the  previous  section.  The  procedure  is  continued  until  it  is  no  longer  possible,  at  which  time  both 
the  index  and  value  of  the  coefficient  which  failed  to  be  condensed  in  this  manner  are  recorded. 
We  iterate  this  procedure  until  it  is  no  longer  efficient.  Our  computational  experiments  with 
images  have  shown  that  a  fixed  value  of  four  iterations  provides  near  optimal  results  on  standard 
test  images,  but  there  are  many  variations  that  may  be  applied  to  tune  this  type  of  algorithm 
depending  on  the  class  of  images  being  analyzed. 

The  interleaving  process  is  simply  to  take  two  signed  integers  (denoted,  for  example,  by  c\  and 
C2)  which  each  require  l  bits  for  their  representation,  and  pack  them  into  21  bits,  but  we  ensure 
that  the  resulting  integer  (which  we  call  cond(c\,  C2))  is  still  relatively  small.  One  method  for 


implementing  this  is  to  use  the  highest  and  lowest  order  bits  of  the  integer  cond(c\ ,  C2)  for  the  sign 
bits  of  ci  and  C2,  respectively.  The  sgn  of  c  is  defined  to  be  ‘O’  if  c  is  nonnegative,  and  ‘T  otherwise. 
Next  we  define  for  any  signed  integer  c,  val(c)  =  |c|  —  sgn(c).  The  positive  integers  val(ci)  and 
val(c 2)  are  packed  into  the  remaining  bits  of  cond{c\,C2)  by  using  the  lowest  odd  bits  for  ci  and 
the  lowest  even  bits  for  02.  It  is  clear  that  if 

-21-1  <  Ci  <  2l~\  (i  =  1,2) 

then 

— 221^1  <  cond(ci,C2 )  <  2ll~1 . 

For  the  purposes  of  illustration,  let  ci  =  2  and  C2  =  —4,  then  sgn(ci)  =  0,  sgn(c2)  =  1  and, 
in  binary  form,  val{c\)  =  010,  val(c2)  =  Oil.  In  this  case,  val  {cond{c\,  C2))  =  00011011  and 
sgn(cond(ci,C2))  =  0. 

This  process  on  roughly-ordered  quantized  coefficients  may  be  summarized  as  follows:  Proceed¬ 
ing  through  the  appropriate  matrices  in  (16),  we  find  the  first  integer  coefficient  which  is  outside 
the  range  [—128,127]  and  record  its  index  ri\ .  We  then  condense  the  pairs  of  coefficients  up  to 
that  index  by  interleaving  adjacent  pairs.  If  rii  is  smaller  than  a  given  value  (i.e.  coarse  enough,  so 
that  few  coefficients  remain),  we  terminate  the  process,  otherwise  we  repeat  the  procedure  on  the 
resulting  data  set  until  stopping  criteria  are  satisfied. 

The  final  step  of  the  algorithm  applies  Q-Coder2  [17]  to  the  encoded  data  as  well  as  to  the 
exceptional  data  and  their  locations,  in  order  to  obtain  a  final  compressed  data  file.  The  recon¬ 
struction  process  is  straightforward  and  well-defined,  since  we  keep  the  record  of  n^.  The  algorithm 
as  described  is  clearly  lossless  up  to  the  original  quantization  procedure.  Variants  of  the  complete 
procedure  also  provide  lossless  compression  of  data  with  progressive  transmission.  The  algorithms 
may  be  further  tuned  using  various  additional  quantization  methods  preserving  the  approximation 
error  and  varying  the  stopping  criteria  according  the  classes  of  image  data. 

5.2  The  Bit-Stream  Method. 

The  “bit-stream”  encoding  method  relies  on  the  spatial  and  frequency  correlations  among  the 
quantized  wavelet  coefficients,  which  are  extracted  by  the  wavelet  decomposition  process,  and 
should  be  expected  to  consistently  perform  better  than  the  interleaving  method  just  described. 
For  simplicity,  we  consider  the  regular  wavelet  basis  and  indicate  how  to  modify  the  algorithm 
for  the  hyperbolic  basis  at  the  end  of  this  section.  We  briefly  describe  the  general  procedure  and 
then  use  an  example  to  illustrate  the  details  of  the  encoding  process.  We  will  begin  with  a  matrix 
A  listing  only  quantized  regular  basis  wavelet  coefficients  (i.e.,  any  of  the  matrices  M5,  M7,  or 
Mg,  corresponding  respectively  to  ipip,  ipip,  or  ip'ip).  From  A  we  generate  a  vector  of  nonzero 
coefficients  c,  and  three  bitstream  linear  arrays  which  encode  the  index  positions:  the  linear  array 
S  will  encode  the  strongly  significant  spatial  correlations  and  implicitly  contains  some  frequency 
correlations,  F  encodes  the  frequency  correlations,  and  R  denotes  the  array  which  encodes  weakly 
significant  spatial  correlations. 

Since  the  majority  of  the  coefficients  of  A  are  zero,  we  write  each  nonzero  coefficient  into  c 
according  to  its  index  position  (ordered  first  by  rows,  next  by  columns,  and  so  on),  after  subtracting 
one  from  each  positive  coefficient  in  order  to  slightly  decrease  its  size  for  efficient  packing.  We  use 

2Our  version  of  Q-coder  was  implemented  by  V.  Zanev  [23]. 
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X  to  denote  the  “indicator  matrix”  of  A,  i.e.  the  characteristic  function  of  the  support  set  of 
indices  of  nonzero  coefficients  of  A.  For  exampie,  if  we  take  the  7x7  matrix: 


0 

0 

1 

-1 

0 

0 

-3 

5 

1 

0 

0 

0 

0 

4 

6 

7 

0 

1 

CO 

2 

1 

5 

the  corresponding  indicator  matrix  X  is  given  by 


o 

0 

1 

1 

0 

0 

1 

1 

1 

0 

0 

0 

o 

1 

1 

1 

0 

1 

1 

1 

1 

and  the  coefficient  matrix  with  positive  coefficients  reduced  by  one  becomes 


0 

0 

0 

-1 

0 

0 

-3 

4 

0 

0 

0 

0 

0 

3 

5 

6 

0 

-3 

1 

0 

4 

which  is  encoded  into  the  array  c,  corresponding  to  the  “nonzero  coefficients”  of  A, 


(23) 


(24) 


(25) 


c  =  (0  —  1  —3  4  0  3  5  6—3  1  0  4). 


Obviousiy,  A  may  be  recovered  from  c  and  the  indicator  matrix  X. 

The  correiations  within  the  matrix  X  are  encoded  by  using  two  types  of  mappings,  one  for 
iocaiiy  encoding  spatiai  correiations  within  a  given  frequency  ievei: 


adjacent  symbol  pairs 
from  matrix  X^' 

encoded  in 
the  matrix  xj*+1^ 

encoded  in  the 

bitstream  array  S 

encoded  in  the 

bitstream  array  R 

0  0 

0 

0  1 

1 

0 

0 

1  0 

1 

0 

1 

1  1 

1 

1 

(26) 


10 


and  another  for  the  correlations  among  adjacent  frequencies: 


symbols  in  submatrices 

encoded  in  the 

encoded  in  the 

(2) 

from  two  levels  of  XT 

matrix  Xy 

array  F 

0  0 

0 

0  1 

1 

0 

1  0 

1 

-1 

1  1 

1 

1 

(27) 


To  illustrate  the  encoding  procedure,  we  continue  with  the  7x7  example  matrix  X  displayed 
in  (24).  We  set  Xo  =  X  and  work  our  way  up  from  the  finest  level  (i.e.  highest  frequency)  to  the 
coarsest  by  a  recursive  procedure.  The  finest  level  indices  for  Xo  are  those  in  its  4x4  upper  left 
submatrix.  We  apply  the  mapping  (26)  first  to  adjacent  pairs  in  the  rows  of  this  submatrix,  to 
obtain 


X?  :=  X0 


0 

0 

1 

1 

0 

0 

1 

1 

1 

0 

0 

0 

0 

1 

1 

1 

0 

1 

1 

1 

1 

o 

1 

o 

1 

1 

0 

— 

1 

0 

1 

1 

1 

1 

S  =  (11001),  R  =  (10). 


We  next  apply  this  procedure  to  the  columns  of  X 
S  and  R  to  obtain 


(i) 

l  > 


appending  to  the  appropriate  bitstreams 


o 

1 

o 

1 

1 

0 

— 

1 

0 

1 

1 

1 

1 

X 


o 

1 

1 

1 

0 

1 

1 

1 

1 

S  =  (11001  110),  R  =  (10  0). 


At  this  stage  we  have  encoded  the  horizontal  and  vertical  spatial  correlations  within  the  current 
frequency  level.  The  next  step  of  the  process  is  to  encode  frequency  information  from  the  current 
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level  to  the  next  coarsest  using  the  mapping  (27).  This  is  accomplished  by  comparing  the  two  2x2 
submatrices,  position  by  position,  from  the  two  levels  (with  finer  first): 


0 

1 

1 

1 

0 

1 

1 

1 

1 

X 


0 

1 

1 

~T 

1 

s 


(11001110),  R  =  (10  0),  F  =  (111). 


( 3) 

The  matrix  Xo  is  now  reduced  (i.e.,  reduced  to  Xi  :=  Xj  ;)  from  having  three  frequency  levels  to 
two,  together  with  the  bitstreams  S,  R,  and  F.  Now  we  iterate  these  three  steps  to  further  reduce 
the  number  of  submatrices  in  X  while  appending  to  the  bitstreams.  It  is  clear  that  the  process  is 
invertible  at  each  stage  and  therefore  the  encoding  is  lossless.  The  full  recursive  process  with  our 
simple  example  is  to  encode  the  horizontal  spatial  correlations  as  step  1, 


0 

1 

1 

~r 

1 

1 

~T 

1 

S  =  (11001  110  01),  R  =  (10  0  0),  F  =  (111), 


to  encode  the  vertical  spatial  correlations  as  step  2, 


1 

~T 

1 

1 

1 

S  =  (11001  110  01  1),  R  =  (10  0  0),  F  =  (111), 


and  to  encode  corresponding  frequency  correlations  in  step  3: 


S  =  (11001  110  01  1), 


R  =  (10  0  0),  F  =  (111  1). 


At  the  last  stage,  we  append  to  S  the  singleton  entry  of  Xn_i  if  n  is  the  number  of  frequency 
levels.  To  reconstruct  the  matrix  A,  we  need  only  the  bitstreams  S,  R,  and  F,  together  with  the 
integer  coefficient  array  c.  The  original  thresholded  matrix  of  wavelet  coefficients  is  obtained  from 
A  and  the  threshold  constant  e.  Standard  arithmetic  encoders  are  used  to  compress  both  c  and 
the  bitstreams  S  and  F.  The  bitstream  R  has  entries  which  are  equally  likely  and  therefore  further 
compression  should  not  be  expected  to  be  effective  on  this  array. 

We  wish  to  emphasize  that  this  algorithm  is  quite  simple  to  extend  in  order  to  handle  higher 

(2) 

dimensional  data.  After  computing  Xj^  which  encodes  horizontal  and  vertical  correlations,  we 
may  apply  the  procedure  again  to  obtain  “depth”  correlations.  Another  simple  modification  is  to 
handle  all  three  matrices  M5,  Mj,  and  Ms,  simultaneously,  by  consolidating  the  scaling-wavelet 
correlations  componentwise  into  an  additional  bitstream  SW.  The  algorithm  is  also  easily  modified 
to  handle  rectangular  data  by  using  a  rectangular  “template”  to  combine  level  to  level  instead  of 
the  square  “template”  described  above.  When  the  dimensions  are  not  a  power  of  two,  we  use  our 
extension  process  to  provide  the  additional  data  to  allow  the  algorithm  to  proceed. 

Finally,  fixed  bases  representations  other  than  “regular”  bases  can  be  handled  with  again  only 
slight  modifications.  The  algorithm  we  have  described  for  the  “regular”  basis  uses  wavelet  (and 


12 


scaling)  basis  elements  of  the  same  frequency  to  determine  spatial  correlations.  Recall  from  section  2 
that  the  “hyperbolic”  basis  consists  of  all  wavelet  coefficients  plus  two  “wings”  from  the  matrices 
M5  and  M-j.  We  consider  at  each  frequency  level  k  those  coefficients  for  which  at  least  one  of  the 
two  frequencies  of  the  tensor  product  factors  is  of  level  k.  We  reduce  this  template  first  in  the 
horizontal  direction  and  then  in  the  vertical  using  the  operation  (26)  as  was  used  for  the  pairs  of 
regular  basis  coefficients.  As  before  we  can  then  compare  the  resulting  reduced  k- th  level  to  next 
coarsest  level  component  by  component  using  the  operation  in  (27). 

We  finish  our  description  of  the  “bit-stream  encoding”  by  stating  a  couple  of  observations 
which  are  not  apparent  in  our  contrived  7x7  example  indicator  matrix  displayed  in  (24).  In  our 
experiments  with  standard  test  images,  the  bitstreams  S  consist  primarily  of  zeros,  since  for  these 
images  the  pair  (11)  appears  infrequently  compared  to  the  pairs  (10)  and  (01)  and  therefore  the 
bitstream  should  compress  significantly.  A  wavelet  coefficient  will  normally  be  large  when  an  edge 
or  sudden  change  occurs.  Depending  on  the  classification  of  the  image,  these  changes  typically 
occur  over  lower  dimensional  sets  which  lead  to  high  occurances  of  the  patterns  (10)  and  (01). 

Shapiro  observed  in  his  quadtree  formulation  for  the  regular  wavelet  basis,  that  typically  if  a 
coefficient  is  zero,  then  with  high  probabilty  the  quadtree  branch  below  that  node  would  consist 
of  zeros.  Our  experiments  have  shown  that  approximately  10%  of  the  time,  if  an  entry  in  the 
coefficient  array  vanishes,  then  the  coefficients  in  the  higher  frequencies  which  would  be  combined 
in  the  reduction  process  with  this  coefficient  (i.e.  whose  indexing  dyadic  squares  are  contained  in 
the  indexing  dyadic  square  of  this  coefficient)  will  all  be  zero.  Since  our  process  essentially  builds 
the  graph  structure  of  the  nonzero  entries,  then  the  zero  entries  are  essentially  ignored,  which 
indicates  that  compression  should  be  effective. 

6.  Computational  Results. 

In  this  section,  we  present  computational  results  comparing  an  implementation  (matching  the 
published  results  in  [20])  of  Sharpiro’s  EZW  algorithm,  the  simplified  interleaving  algorithm,  and 
the  bit-stream  encoding  algorithm  described  in  this  paper.  We  apply  these  algorithms  to  both 
regular  and  hyperbolic  wavelets  using  the  6-10  biorthogonal  filters  (see  Table  I  for  the  coefficients). 
In  each  of  these  experiments,  the  resulting  files  are  passed  through  the  Q-coder  algorithm  to  finalize 
the  compression.  We  have  chosen  four  standard  test  images  (see  Figure  1)  in  order  to  compare 
the  compression  rates  and  the  PSNR  (Peak  Signal  to  Noise  Ratio)  of  the  methods.  The  PSNR  is 
considered  to  be  a  measure  of  image  quality  and  is  given  by  the  formula 

(255  \ 

- — - 

II/-/IIJ 

where  /  is  the  original  image,  and  /  is  the  reconstructed  image  from  the  compressed  data,  and 
£2  is  the  least  squares  norm.  The  four  test  images  range  from  relatively  smooth  to  rough  images 
(Lena,  Barbara,  City,  Destroyer).  All  experiments  are  performed  by  encoding  and  then  decoding 
512  x  512  square  images  for  comparison.  As  mentioned  above,  we  have  implemented  our  own  version 
of  Shapiro’s  EZW  algorithm  which  performs  according  to  the  published  results  in  [20],  while  the 
general  Q-coder  was  implemented  by  our  colleague  Zanev  [23]  as  described  in  [17]. 

We  argue  that  our  algorithms,  using  both  regular  and  hyperbolic  wavelets,  are  at  least  com¬ 
parable  to  Shapiro’s  EZW  algorithm  and  are  easy  to  implement  in  both  software  and  hardware. 
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Furthermore,  the  methods  are  efficient  and  easily  extendible  to  arbitrary  m  x  n  rectangular  images 
index  limits  and,  perhaps  more  importantly,  to  higher  dimensions. 

In  Table  II,  we  compare  compression  rates  of  the  Lena  test  image  for  Shapiro’s  method  in 
[20],  our  implementation  of  his  method,  our  elementary  interleaving  encoding,  and  our  bit-stream 
encoding  using  both  regular  and  hyperbolic  wavelets.  Table  III  provides  similar  information  about 
the  Barbara  test  image.  These  two  images  were  selected  since  they  are  standard  and  were  used 
as  test  images  in  [20].  The  reasoning  behind  the  inclusion  of  the  results  from  our  implementation 
of  Shapiro’s  method  is  to  show  that  it  performs  consistent  with  published  results  in  [20],  for  both 
the  Lena  and  Barbara  images.  We  can  then  get  an  indication  of  the  performance  of  that  method 
on  other  less  smooth  test  images  for  which  test  results  are  not  available  in  [20].  Tables  IV- V 
compares  our  implementation  of  Shapiro’s  method  against  the  encoding  applied  to  both  Regular 
and  Hyperbolic  biorthogonal  wavelet  decompositions  for  the  City  and  Destroyer  test  images.  The 
computational  results  for  bit-stream  encoding  against  Shapiro’s  method  is  also  presented  in  Plots  1— 
3.  The  results  of  Table  II  are  presented  in  Plot  1  which  plots  the  PSNR  against  compression  rate, 
which  ranges  from  8  to  256.  Plots  2-3  presents  the  corresponding  information  for  the  Barbara, 
City  and  Destroyer  test  images,  respectively.  It  is  clear  from  these  plots  and  tables  that  all  six 
methods  are  comparable  with  some  advantage  of  the  interleaving  encoding  for  the  Destroyer  image. 
Figure  2  compares  the  compressed  Lena  images  by  the  five  methods  against  the  original  image  using 
a  compression  rate  of  128.  Figure  3  performs  the  same  role  for  the  rougher  Barbara  image.  These 
figures  give  a  qualitative  measure  of  compression  for  the  three  methods  holding  the  compression 
ratio  fixed. 

The  images  we  use  were  also  selected  since  they  have  varying  degrees  of  smoothness  in  Besov 
spaces  and  the  nonlinear  approximation  theory  of  DeVore  [10]  provides  a  basis  for  analysis  of  these 
results.  From  the  tables  and  plots  we  can  see  that  generally  the  bit-stream  method  gives  the 
highest  compression  while  the  interleaving  method  gives  the  lowest,  with  the  exception  for  the 
Destroyer  image.  It  is  well  known  that  the  regular  basis  works  better  than  the  hyperbolic  basis 
for  the  ‘smooth’  Lena  image,  while  hyperbolic  basis  is  a  better  choice  for  ‘rougher’  Barbara  image. 
In  combination  with  our  algorithms,  both  interleaving  and  bit-stream  encoding,  this  property  is 
retained. 

To  indicate  the  strong  potential  of  using  the  encoders  and  multidimensional  biorthogonal 
wavelets  for  data  compression,  we  apply  these  methods  to  a  relatively  rough  contaminant  con¬ 
centration  front  invading  the  groundwater  in  a  porous  media  (see  Figure  4).  The  scalar  data  field 
is  evaluated  on  a  nonuniform  grid  with  137  x  129  x  11  vertices.  In  this  experiment,  we  use  the 
hyperbolic  basis  with  the  symmetric  6-10  filter  in  each  coordinate  and  the  interleaving  method  for 
encoding  the  significant  coefficients.  We  have  rendered  this  scalar  concentration  field  according  to 
the  color  map  located  at  the  bottom  of  each  image.  In  this  figure,  we  have  compared  the  renderings 
of  the  original  file  against  those  reconstructed  from  compressed  files  with  compression  rates  of  172, 
493,  and  1200,  respectively.  Similar  compression  rates  have  been  obtained  for  compression  of  the 
grid  data,  with  no  visual  loss  of  information.  In  practice,  the  original  grid  data  are  generated 
from  scattered  data  whose  measurements  contain  varying  degrees  of  error.  One  may  then  perhaps 
argue  that  the  compressed  (denoised)  grid  provides  at  least  as  good  a  representation  of  the  physical 
grid  as  the  original  file.  These  same  algorithms  have  been  applied  to  video  and  to  time-dependent 
three-dimensional  fields  and  will  be  reported  on  elsewhere. 
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7.  Remarks 


As  stated  in  the  introduction,  our  original  intent  for  this  study  was  the  development  of  reli¬ 
able  data  compression  algorithms  for  three-dimensional  scalar  and  vector  fields  with  general  index 
limits  (i.e.  non-dyadic)  for  use  in  interactive  steering  of  simulations  on  remote  massively  parallel 
supercomputers  [14,  15].  In  this  report  we  have  described  two  methods  for  encoding  the  biorthog- 
onal  wavelet  coefficients  and  have  demonstrated  that  they  generally  outperform  other  compression 
methods.  We  also  observe  that  bit-stream  encoding  is  typically  the  superior  of  the  two  methods 
for  most  images.  Our  algorithms,  though  primarily  tested  for  2  dimensional  square  images,  have 
natural  extensions  to  higher  dimensions  for  any  logically  rectangular  data.  These  algorithms  also 
are  well  suited  for  use  with  adaptive  bases,  for  example  the  minimal  entropy  basis  [4]  and  adaptive 
nonlinear  bases  [18]. 

Shannon’s  entropy  may  be  used  to  understand  the  effectiveness  of  the  bitstream  encoding.  If  S  is 
a  random  variable  which  takes  k  values  vi,  v2, ■« . . ,  with  corresponding  probabilities  pi,p2,  ■  ■  ■  ,Pk, 
then  the  entropy  of  S  may  be  defined  as 

k 

H{S)  =  -J2  pi  log 2pi 
i= 1 

and  is  considered  a  measure  of  the  complexity  of  S.  The  Noiseless  Coding  Theorem  of  information 
theory  [1,  19]  guarantees  that  the  minimal  expected  length  of  encoded  data  among  all  uniquely 
decipherable  binary  encoding  schemes  may  essentially  be  obtained  by  an  extension  of  Huffman 
encoding.  Shannon’s  theory  states  that  if  the  size  of  the  data  is  n,  then  nH(S)  is  the  expected 
number  of  bits  required  to  represent  S  and  is  then  called  the  complexity  of  S.  This  branch  of 
information  theory  also  provides  prefix  encoders,  such  as  Huffman  encoding,  which  on  average 
provide  minimal  length  encoding  of  the  data  (i.e.  nH(S)).  Datum  in  the  data  sets  are  treated 
as  values  generated  by  repeated  application  of  a  random  variable,  that  is,  it  is  assumed  that  the 
values  are  produced  by  independent,  identically  distributed  random  variables.  In  practice,  the 
probabilities  are  then  estimated  by  the  frequency  table  of  the  data.  The  bitstream  preprocessor 
may  be  regarded  as  a  prefix  encoder  which  extracts  and  encodes  significant  wavelet  correlations  by 
means  of  a  local  space-frequency  Huffman  procedure  (see  definitions  (26)  and  (27))  and  thereby 
prepares  the  data  for  an  enhanced  application  of  standard  data  compression  algorithms  such  as 
that  which  is  incorporated  in  Q-coder. 

Theoretically,  we  may  describe  the  performance  of  the  bitstream  algorithm  as  follows.  Given  a 
binary  data  set  S  with  complexity  nH(S),  if  we  apply  one  step  bitstream  encoding  on  S ,  then  S 
will  be  represented  as  a  4-ary  data  set  S'  whose  size  is  §.  The  complexity  of  S'  is  Under 

these  conditions,  one  can  show  [13]  that  <  nH(S)  and  that  bitstream  encoding  decreases 

the  complexity  of  the  original  data.  This  result  not  only  assures  the  effectiveness  of  the  bitstream 
encoding,  but  also  provides  a  parameter  [13],  which  depends  on  relative  ratios  of  the  frequencies  of 
the  four  symbols,  to  measure  its  performance. 

Typical  data  which  we  have  encountered  in  large  scale  scientific  simulations  and  medical  im¬ 
age  processing  may  be  quite  sizable  and  has  motivated  further  study  of  parallel  versions  of  our 
algorithms,  as  well  as  progressive  transmission.  These  efforts  will  be  reported  on  elsewhere. 
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1 

a.  I 

A 

-4 

0 

0.01337437 

-3 

0 

0.00494231 

-2 

-0.09127176 

-0.04754360 

-1 

0.03372823 

0.09432042 

0 

0.55754352 

0.43490656 

1 

0.55754352 

0.43490656 

2 

0.03372823 

0.09432042 

3 

-0.09127176 

-0.04754360 

4 

0 

0.00494231 

5 

0 

0.01337437 

TABLE  I.  6-10  Biorthogonal  Wavelet  Filter  Coefficients  for  Com¬ 
putational  Experiments  (see  equation  (13)). 


Compression 

Rate 

Shapiro 

PSNR  (dB) 

use 

Shapiro 
PSNR  (dB) 

Regular 
Condensation 
PSNR  (dB) 

Hyperbolic 
Condensation 
PSNR  (dB) 

Regular 

Bit  Stream 
PSNR  (dB) 

Hyperbolic 
Bit  Stream 
PSNR  (dB) 

8 

39.55 

39.16 

39.22 

38.93 

39.57 

39.24 

16 

36.28 

36.05 

35.88 

35.82 

36.65 

36.33 

32 

33.17 

32.83 

33.15 

32.82 

33.72 

33.45 

64 

30.23 

30.01 

29.91 

29.78 

30.89 

30.51 

128 

27.54 

27.66 

27.66 

27.16 

28.21 

27.90 

256 

25.38 

25.33 

25.26 

25.07 

25.86 

25.53 

TABLE  II.  Quantitative  Results  Comparing  Peak-Signal  to  Noise  Ra¬ 
tios  (PSNR)  for  512  x  512  Lena  Image  Using  Various  Meth¬ 
ods  with  Compression  Rates  from  8-256. 


Compression 

Rate 

Shapiro 

PSNR  (dB) 

use 

Shapiro 
PSNR  (dB) 

Regular 
Condensation 
PSNR  (dB) 

Hyperbolic 
Condensation 
PSNR  (dB) 

Regular 

Bit  Stream 
PSNR  (dB) 

Hyperbolic 
Bit  Stream 
PSNR  (dB) 

8 

35.14 

33.64 

35.10 

35.20 

35.59 

35.65 

16 

30.53 

29.09 

30.55 

30.49 

30.81 

30.96 

32 

26.77 

26.09 

26.66 

27.04 

27.22 

27.50 

64 

24.03 

24.22 

24.05 

24.18 

24.43 

24.69 

128 

23.10 

22.76 

22.83 

22.71 

23.10 

23.00 

256 

21.94 

21.92 

21.68 

21.50 

22.14 

21.92 

TABLE  III.  Quantitative  Results  Comparing  Peak-Signal  to  Noise  Ra¬ 
tios  (PSNR)  for  512  x  512  Barbara  Image  Using  Various 
Methods  with  Compression  Rates  from  8-256. 


Compression 

Rate 

use 

Shapiro 
PSNR  (dB) 

Regular 
Condensation 
PSNR  (dB) 

Hyperbolic 
Condensation 
PSNR  (dB) 

Regular 

Bit  Stream 
PSNR  (dB) 

Hyperbolic 
Bit  Stream 
PSNR  (dB) 

8 

24.04 

24.72 

24.19 

25.02 

24.55 

16 

22.03 

22.25 

21.81 

22.43 

21.98 

32 

20.33 

20.63 

20.26 

20.77 

20.39 

64 

19.14 

19.42 

19.17 

19.53 

19.24 

128 

18.40 

18.44 

18.26 

18.57 

18.40 

256 

17.60 

17.70 

17.50 

17.81 

17.61 

TABLE  IV.  Quantitative  Results  Comparing  Peak-Signal  to  Noise  Ra¬ 
tios  (PSNR)  for  512  x  512  City  Image  Using  Various  Meth¬ 
ods  with  Compression  Rates  from  8-256. 


Compression 

Rate 

use 

Shapiro 
PSNR  (dB) 

Regular 
Condensation 
PSNR  (dB) 

Hyperbolic 
Condensation 
PSNR  (dB) 

Regular 

Bit  Stream 
PSNR  (dB) 

Hyperbolic 
Bit  Stream 
PSNR  (dB) 

8 

36.89 

40.85 

40.67 

40.38 

40.56 

16 

33.31 

34.78 

35.35 

34.62 

35.31 

32 

31.42 

32.08 

32.50 

31.97 

32.51 

64 

29.15 

30.69 

30.86 

30.78 

30.92 

128 

26.26 

28.39 

28.40 

28.48 

28.52 

256 

24.10 

25.64 

25.88 

25.70 

26.23 

TABLE  V.  Quantitative  Results  Comparing  Peak-Signal  to  Noise  Ra¬ 
tios  (PSNR)  for  512  x  512  Destroyer  Image  Using  Various 
Methods  with  Compression  Rates  from  8-256. 


Figure  1.  Original  test  images  used  in  com] 
(c)  City,  and  (d)  Destroyer. 


experiments:  (a)  Lena,  (b)  Barbara, 
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Figure  3.  Compression  of  Barbra  Image  with  a  compression  rate  of  128  comparing  methods:  (a) 
original  image,  (b)  Shapiro  (South  Carolina  implementation),  (c)  Hyperbolic  wavelets 
with  interleaving,  (d)  Regular  wavelets  with  interleaving,  (e)  Hyperbolic  wavelets  with 
bit- stream  and  (f)  Regular  wavelets  with  bit- stream. 


(c) 


(d) 


Figure  4.  Comparision  of  3D  data  compression,  using  the  interleaving  encoding,  to  orig¬ 
inal  3D  scalar  field  with  varying  compression  rates:  (a)  Original  data,  (b)  172 
times  compression,  (c)  493  times  compression,  and  (d)  1200  times  compression. 
All  images  are  rendered  with  the  same  orthogonal  slices  and  an  isosurface  of 
.3  of  maximum  concentration  percentage. 


Plot  1.  Plots  of  PS  NR  (Peak  Signal  to  Noise  Ratio)  as  a  function  of  the  compression 
rate  for  the  Lena  test  image  using  Shapiro’s  EZM  algorithm  (dotted),  the  USC 
version  of  Shapiro  EZM  algorithm  (solid)  and  the  elementary  encoding  applied 
to  regular  (dashed)  and  hyperbolic  (dash-dotted)  wavelets. 


Plot  2.  Plots  of  PS  NR  (Peak  Signal  to  Noise  Ratio)  as  a  function  of  the  compression 
rate  for  the  Barbara  test  image  using  Shapiro’s  EZM  algorithm  (dotted),  the 
USC  version  of  Shapiro  EZM  algorithm  (solid)  and  the  elementary  encoding 
applied  to  regular  (dashed)  and  hyperbolic  (dash-dotted)  wavelets. 


Compression  Rate 


Compression  Rate 


(a) 


(b) 


Plot  3.  Plots  of  PSNR  (Peak  Signal  to  Noise  Ratio)  as  a  function  of  the  compression 
rate  using  the  USC  version  of  Shapiro’s  method  (solid)  and  our  encoding  algo¬ 
rithm  applied  to  regular  (dashed)  and  hyperbolic  (dash-dotted)  wavelets  for: 
(a)  City  image  and  (b)  Destroyer  image. 


